package scu.maqiang.fes;

import static scu.maqiang.numeric.Constants.GP3;
import static scu.maqiang.numeric.Constants.GW3;

import java.util.Arrays;

import scu.maqiang.mesh.Mesh;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.DoubleIntegerPair;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.Pair;
import scu.maqiang.numeric.ParamCheck;
import scu.maqiang.numeric.SRMatrix;

public class FES2Q81 extends FES {
	
	public FES2Q81(Mesh mesh2D) {
		super(mesh2D, 1);
		eleStiff = FES2Q81::elementStiff;
		eleMass = FES2Q81::elementMass;
		eleSource = FES2Q81::elementSource;
	}
	
	public double getValue(double[] u, double[] xyz) {
		ParamCheck.checkEqual(u.length, this.ndof);
		int[] eleIdx = new int[nPerEle];
		double[] eleVal = new double[nPerEle];
	    double[] uu = new double[8];
	    double value = 0;
	    double[][] coord = new double[8][2];
	    for(int i = 0, n = mesh.getNt(); i < n; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	if (CG.pointInQ4(coord, xyz)) {
	    		CG.Q8ShapeFunction(coord, xyz, eleVal);
	    		for(int j = 0; j < 8; j++) {
	    			uu[j] = u[ele[j]];
	    			eleIdx[j]= ele[j];
	    			value += uu[j] * eleVal[j];
	    		}
	    		return value;
	    	}	    	
	    }
	    throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
	}
	
	public double getValue(double[] u, double[] xyz, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
	    double[] uu = new double[8];
	    double value = 0;
	    double[][] coord = new double[8][2];
	    for(int i = 0, n = mesh.getNt(); i < n; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	if (CG.pointInQ4(coord, xyz)) {
	    		CG.Q8ShapeFunction(coord, xyz, eleVal);
	    		for(int j = 0; j < 8; j++) {
	    			uu[j] = u[ele[j]];
	    			eleIdx[j]= ele[j];
	    			value += uu[j] * eleVal[j];
	    		}
	    		return value;
	    	}	    	
	    }
	    throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
	}
	
	public Pair<Double, Integer> getValue(double[] u, double[] xy, int previousEleIdx, double[][] coef, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double value = 0;
		int[] ele = mesh.getElement(previousEleIdx);
		double x = xy[0], y = xy[1];
		double[] v = new double[] {1.0, x, y, x * x, x * y, y * y, x * x * y, x * y * y};
		mesh.getCoordPerElement(ele, coordRegion);
		if (CG.pointInQ4(coordRegion, xy)) {
			MVO.matmul(coef, v, eleVal);
			for (int j = 0; j < 8; j++) {
				value += u[ele[j]] * eleVal[j];
			}
			return new Pair<Double, Integer>(value, previousEleIdx);
		} else {
			for (int i = 0, n = mesh.getNt(); i < n; i++) {
				ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				if (CG.pointInQ4(coordRegion, xy)) {
					//CG.Q8ShapeFunction(coord, xy, eleVal);
					CG.Q8ShapeFunctionCoef(coordRegion, coef);
					MVO.matmul(coef, v, eleVal);
					for (int j = 0; j < 8; j++) {
						eleIdx[j] = ele[j];
						value += u[ele[j]] * eleVal[j];
					}
					return new Pair<Double, Integer>(value, i);
				}
			}
			throw new IllegalArgumentException("The Point ( " + Arrays.toString(xy) + " ) is not in the domain!");
		}

	}
	
	public void getValue(double[] u, double[] xy, DoubleIntegerPair data, double[][] coef, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double value = 0;
		int[] ele = mesh.getElement(data.second);
		double x = xy[0], y = xy[1];
		double[] v = new double[] {1.0, x, y, x * x, x * y, y * y, x * x * y, x * y * y};
		mesh.getCoordPerElement(ele, coordRegion);
		if (CG.pointInQ4(coordRegion, xy)) {
			MVO.matmul(coef, v, eleVal);
			for (int j = 0; j < 8; j++) {
				value += u[ele[j]] * eleVal[j];
			}
			//data.setFirst(value);
			data.first = value;
			return;
			//return new Pair<Double, Integer>(value, previousEleIdx);
		} else {
			for (int i = 0, n = mesh.getNt(); i < n; i++) {
				ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				if (CG.pointInQ4(coordRegion, xy)) {
					//CG.Q8ShapeFunction(coord, xy, eleVal);
					CG.Q8ShapeFunctionCoef(coordRegion, coef);
					MVO.matmul(coef, v, eleVal);
					for (int j = 0; j < 8; j++) {
						eleIdx[j] = ele[j];
						value += u[ele[j]] * eleVal[j];
					}
					//data.setFirst(value);
					//data.setSecond(i);
					data.first = value;
					data.second = i;
					return;
					//return new Pair<Double, Integer>(value, i);
				}
			}
			throw new IllegalArgumentException("The Point ( " + Arrays.toString(xy) + " ) is not in the domain!");
		}

	}
	
	public void getValue_Square(double[] u, double[] xy, DoubleIntegerPair data, double[][] coef, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double value = 0;
		int[] ele = mesh.getElement(data.second);
		double x = xy[0], y = xy[1];
		double[] v = new double[] {1.0, x, y, x * x, x * y, y * y, x * x * y, x * y * y};
		mesh.getCoordPerElement(ele, coordRegion);
		if (CG.pointInSquare(coordRegion, xy)) {
			MVO.matmul(coef, v, eleVal);
			for (int j = 0; j < 8; j++) {
				value += u[ele[j]] * eleVal[j];
			}
			//data.setFirst(value);
			data.first = value;
			return;
			//return new Pair<Double, Integer>(value, previousEleIdx);
		} else {
			for (int i = 0, n = mesh.getNt(); i < n; i++) {
				ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coordRegion);
				if (CG.pointInSquare(coordRegion, xy)) {
					//CG.Q8ShapeFunction(coord, xy, eleVal);
					CG.Q8ShapeFunctionCoef(coordRegion, coef);
					MVO.matmul(coef, v, eleVal);
					for (int j = 0; j < 8; j++) {
						eleIdx[j] = ele[j];
						value += u[ele[j]] * eleVal[j];
					}
					//data.setFirst(value);
					//data.setSecond(i);
					data.first = value;
					data.second = i;
					return;
					//return new Pair<Double, Integer>(value, i);
				}
			}
			throw new IllegalArgumentException("The Point ( " + Arrays.toString(xy) + " ) is not in the domain!");
		}

	}

	/**
	 * ���㵥Ԫ�նȾ���
	 * 
	 * @param xy
	 * @param kxx
	 * @param tp
	 * @param eleMatrix
	 */
	public static void elementStiff(int[] ele, double[][] xy, double[][] kxx, BVPType tp, double[][] eleMatrix) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] v = new double[8];
		double xi, et, detJ, wxarea, c;
		MVO.fill(eleMatrix, 0.0);
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					c = MVO.dot_product(kxx[0], v);
					wxarea = GW3[i] * GW3[j] * detJ * c;
					for (int si = 0; si < 8; si++) {
						for (int sj = 0; sj < 8; sj++) {
							eleMatrix[si][sj] += (dx[si] * dx[sj] + dy[si] * dy[sj]) * wxarea;
						}
					}
				}
			}
			break;
		case AXISYMMETRIC:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					c = MVO.dot_product(kxx[0], v) * MVO.dot_product(x, v);
					wxarea = GW3[i] * GW3[j] * detJ * c;
					for (int si = 0; si < 8; si++) {
						for (int sj = 0; sj < 8; sj++) {
							eleMatrix[si][sj] += (dx[si] * dx[sj] + dy[si] * dy[sj]) * wxarea;
						}
					}
				}
			}
			break;
		default:
			break;
		}
		//System.out.println(MVO.toString(eleMatrix));
	}

	public static void elementMass(int[] ele, double[][] xy, double[][] coef, BVPType tp, double[][] eleMatrix) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] v = new double[8];
		double xi, et, detJ, wxarea, c;
		MVO.fill(eleMatrix, 0.0);
		switch (tp) {
		case CONSIST_MASS:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					c = MVO.dot_product(coef[0], v);
					wxarea = GW3[i] * GW3[j] * detJ * c;
					for (int si = 0; si < 8; si++) {
						for (int sj = 0; sj < 8; sj++) {
							eleMatrix[si][sj] += v[si] * v[sj] * wxarea;
						}
					}
				}
			}
			break;
		case LUMP_MASS:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					c = MVO.dot_product(coef[0], v);
					wxarea = GW3[i] * GW3[j] * detJ * c;
					for (int si = 0; si < 8; si++) {
						eleMatrix[si][si] += v[si] * wxarea;
					}
				}
			}
			break;
		case AXIS_CONSIST_MASS:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					c = MVO.dot_product(coef[0], v) * MVO.dot_product(x, v);
					wxarea = GW3[i] * GW3[j] * detJ * c;
					for (int si = 0; si < 8; si++) {
						for (int sj = 0; sj < 8; sj++) {
							eleMatrix[si][sj] += v[si] * v[sj] * wxarea;
						}
					}
				}
			}
			break;
		case AXIS_LUMP_MASS:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					c = MVO.dot_product(coef[0], v) * MVO.dot_product(x, v);
					wxarea = GW3[i] * GW3[j] * detJ * c;
					for (int si = 0; si < 8; si++) {
						eleMatrix[si][si] += v[si] * wxarea;
					}
				}
			}
			break;
		default:
			break;
		}
	}

	
	public void assembleHeatThetaScheme(ScalarFunc[] coefFunc, double[] param, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
		double[][] kxx = new double[1][nPerEle];
		double[][] rho = new double[1][nPerEle];
		double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
		double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
		double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
		double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
		double dt = param[0];
		double theta = param[1];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			int label = mesh.getElementLabel(i);
			for (int j = 0, nj = nPerEle; j < nj; j++) {
				kxx[0][j] = coefFunc[0].action(coordRegion[j], label, param);
				rho[0][j] = coefFunc[1].action(coordRegion[j], label, param);
			}
			elementStiff(ele, coordRegion, kxx, tp[0], eleMatrixStiff);
			elementMass(ele, coordRegion, rho, tp[1], eleMatrixMass);
			for(int ii = 0; ii < dofPerElement; ii++ ) {
				for(int jj = 0; jj < dofPerElement; jj++) {
					eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
					eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
				}
			}
			LMatrix.assemble(ele, eleLMatrix);
			RMatrix.assemble(ele, eleRMatrix);
		}
	}
	
	public void assembleHeatThetaScheme(double[] coefFunc, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
		double[][] kxx = new double[1][nPerEle];
		double[][] rho = new double[1][nPerEle];
		double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
		double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
		double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
		double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
		Arrays.fill(kxx[0], coefFunc[0]);
		Arrays.fill(rho[0], coefFunc[1]);
		double dt = coefFunc[2];
		double theta = coefFunc[3];		
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			elementStiff(ele, coordRegion, kxx, tp[0], eleMatrixStiff);
			elementMass(ele, coordRegion, rho, tp[1], eleMatrixMass);
//			System.out.println(MVO.toString(eleMatrixMass));
//			System.out.println(MVO.toString(eleMatrixStiff));
//			System.exit(0);
			for(int ii = 0; ii < dofPerElement; ii++ ) {
				for(int jj = 0; jj < dofPerElement; jj++) {
					eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
					eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
				}
			}
			LMatrix.assemble(ele, eleLMatrix);
			RMatrix.assemble(ele, eleRMatrix);
		}
	}

	public static void elementSource(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] v = new double[8];
		double xi, et, detJ, wxarea, c;
		MVO.setValue(eleF, 0.0);
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);

					c = MVO.dot_product(f[0], v);
					wxarea = GW3[i] * GW3[j] * detJ * c;
					for (int si = 0; si < 8; si++) {
						eleF[si] += v[si] * wxarea;
					}
				}
			}
			break;
		case AXISYMMETRIC:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
					c = MVO.dot_product(f[0], v) * MVO.dot_product(v, x);
					wxarea = GW3[i] * GW3[j] * detJ * c;
					for (int si = 0; si < 8; si++) {
						eleF[si] += v[si] * wxarea;
					}
				}
			}
			break;
		default:
			break;
		}
	}
	
	public void computeGradient2(double[] u, double[][] gradU) {
		int len = u.length;
		ParamCheck.checkEqual(len, this.ndof);
		ParamCheck.checkEqual(len, gradU[0].length);
		ParamCheck.checkEqual(len, gradU[1].length);
		
		double[] angle = new double[8];
		//double[][] coord = new double[8][2];
		double[][] grad = new double[8][2];
		double[] angleCount = new double[this.GetNdof()];
		double[] eleU = new double[8];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			for(int j = 0; j < 8; j++) {
				eleU[j] = u[ele[j]];
			}
			CG.gradientValue2DQ8(coordRegion, eleU, grad);
			CG.angle2DQ8(coordRegion, angle);
			for(int j = 0; j < 8; j++) {
				gradU[0][ele[j]] += grad[j][0] * angle[j];
				gradU[1][ele[j]] += grad[j][1] * angle[j];
				angleCount[ele[j]] += angle[j];
			}
		}
		for(int i = 0, n = mesh.getNv(); i < n; i++) {
			gradU[0][i] /= angleCount[i];
			gradU[1][i] /= angleCount[i];
		}
	}
	
	public void computeSecondDerivative(double[] u, double[][] sdU) {
		int len = u.length;
		ParamCheck.checkEqual(len, this.ndof);
		ParamCheck.checkEqual(len, sdU[0].length);
		ParamCheck.checkEqual(len, sdU[1].length);
		ParamCheck.checkEqual(len, sdU[2].length);
		
		double[] angle = new double[8];
		//double[][] coord = new double[8][2];
		double[][] eleSDU = new double[8][3];
		double[] angleCount = new double[this.GetNdof()];
		double[] eleU = new double[8];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			for(int j = 0; j < 8; j++) {
				eleU[j] = u[ele[j]];
			}
			CG.secondDerivativeValue2DQ8(coordRegion, eleU, eleSDU);
			CG.angle2DQ8(coordRegion, angle);
			for(int j = 0; j < 8; j++) {
				sdU[0][ele[j]] += eleSDU[j][0] * angle[j];
				sdU[1][ele[j]] += eleSDU[j][1] * angle[j];
				sdU[2][ele[j]] += eleSDU[j][2] * angle[j];
				angleCount[ele[j]] += angle[j];
			}
		}
		for(int i = 0, n = mesh.getNv(); i < n; i++) {
			sdU[0][i] /= angleCount[i];
			sdU[1][i] /= angleCount[i];
			sdU[2][i] /= angleCount[i];
		}
	}

	public void computeFirstAndSecondDerivative(double[] u, double[][] fdU, double[][] sdU) {
		int len = u.length;
		ParamCheck.checkEqual(len, this.ndof);
		ParamCheck.checkEqual(len, sdU[0].length);
		ParamCheck.checkEqual(len, sdU[1].length);
		ParamCheck.checkEqual(len, sdU[2].length);
		
		double[] angle = new double[8];
		//double[][] coord = new double[8][2];
		double[][] eleFDU = new double[8][2];
		double[][] eleSDU = new double[8][3];
		double[] angleCount = new double[this.GetNdof()];
		double[] eleU = new double[8];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			for(int j = 0; j < 8; j++) {
				eleU[j] = u[ele[j]];
			}
			CG.firstAndSecondDerivativeValue2DQ8(coordRegion, eleU, eleFDU, eleSDU);
			CG.angle2DQ8(coordRegion, angle);
			for(int j = 0; j < 8; j++) {
				fdU[0][ele[j]] += eleFDU[j][0] * angle[j];
				fdU[1][ele[j]] += eleFDU[j][1] * angle[j];
				sdU[0][ele[j]] += eleSDU[j][0] * angle[j];
				sdU[1][ele[j]] += eleSDU[j][1] * angle[j];
				sdU[2][ele[j]] += eleSDU[j][2] * angle[j];
				angleCount[ele[j]] += angle[j];
			}
		}
		for(int i = 0, n = mesh.getNv(); i < n; i++) {
			fdU[0][i] /= angleCount[i];
			fdU[1][i] /= angleCount[i];
			sdU[0][i] /= angleCount[i];
			sdU[1][i] /= angleCount[i];
			sdU[2][i] /= angleCount[i];
		}
	}
}
